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ABSTRACT 

Observations of structure in circumstellar debris discs provide circumstantial evidence 
for the presence of massive planets at large (several tens of AU) orbital radii, where 
the timescale for planet formation via core accretion is prohibitively long. Here, we 
investigate whether a population of distant planets can be produced via outward 
migration subsequent to formation in the inner disc. Two possibilities for significant 
outward migration are identified. First, cores that form early at radii a ~ 10 AU 
can be carried to larger radii via gravitational interaction with the gaseous disc. This 
process is efficient if there is strong mass loss from the disc - either within a cluster 
or due to photoevaporation from a star more massive than the Sun - but does not 
require the extremely destructive environment found, for example, in the core of the 
Orion Nebula. We find that, depending upon the disc model, gas disc migration can 
yield massive planets (several Jupiter masses) at radii of around 20-50 AU. Second, 
interactions within multiple planet systems can drive the outer planet into a large, 
normally highly eccentric orbit. A series of scattering experiments suggests that this 
process is most efficient for lower mass planets within systems of unequal mass ratio. 
This mechanism is a good candidate for explaining the origin of relatively low mass 
giant planets in eccentric orbits at large radii. 

Key words: accretion, accretion discs — stars: formation — stars: pre-main-sequence 
— planetary systems: protoplanetary discs — planets and satellites: formation 
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With one exception (Konacki et al. 2003), all confirmed ex- 
trasolar planets have been discovered by the Doppler veloc- 
ity technique. The selection effects inherent to radial velocity 
surveys (Cumming, Marcy & Butler 1999) favor the detec- 
tion of planets at small orbital radii. To date, about half of 
the known planets have semi-major axis a < 1 AU, while 
the most distant - 55 Cnc d - lies at 5.9 AU from its parent 
star. . Indirect evidence, however, suggests that there could 
be a sizable population of massive planets at much greater 
radii. Recent observations of dusty debris around Vega have 
been interpreted as suggesting the presence of a planet of a 
few Jupiter masses with a > 30 (Wilner et al. 2002). Fur- 
ther, simulations modeling circumstellar dust discs suggest 
a planet lies at a distance of 55 — 65 AU from Epsilon Eridani 
(Ozernoy et al. 2000). 

Forming planets in these outer locations is difficult. Gas 
giants must form before the disc dissipates, at timescales 
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1 From the online Extrasolar Planets Encyclopedia, at 
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no greater than 5 — 10 Myr (Haisch, Lada, & Lada 2001). 
In standard core accretion models (Safranov 1969), the 
timescale for building the core of a giant planet increases 
rapidly with radius, with a tj orm scaling approximately as 
o 2 (Pollack et al. 1996). Although such models are undoubt- 
edly oversimplified (Pollack et al. 1996; Bryden, Lin & Ida 
2000), it is hard to avoid the conclusion that forming mas- 
sive planets at radii of several tens of AU within 10 Myr is 
difficult. Indeed, this has led to the suggestion that Uranus 
and Neptune may have formed at smaller radii in our own 
Solar System (Thommes, Duncan & Levison 1999, 2002). 
Motivated by these issues, we investigate the possibility of 
forming massive planets at small a, followed by outward mi- 
gration. In Sections 2 and 3, we consider sequentially the two 
mechanisms that have been extensively studied in the con- 
text of inward migration: planet-disc interactions (Goldreich 
& Tremaine 1980; Lin & Papaloizou 1986; Lin, Bodenheimer 
& Richardson 1996; Trilling et al. 1998) and gravitational 
scattering after disc dissipation (Rasio & Ford 1996; Weiden- 
schilling & Marzari 1996; Lin & Ida 1997; Ford, Havlickova & 
Rasio 2001; Terquem & Papaloizou 2002). Our conclusions 
are briefly summarized in Section 4. 
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2 MIGRATION VIA GAS DISC 
INTERACTIONS 

2.1 Methods 

We calculate the orbital evolution of massive planets embed- 
ded within an evolving protoplanetary disc using a variant of 
the approach described by Armitage et al. (2002). We use a 
simple, one-dimensional (i.e. vertically averaged) treatment 
to model the evolution of a protoplanetary disc evolving un- 
der the action of both internal viscous torques and external 
torques from one or more embedded planets (Goldreich & 
Tremaine 1980; Lin & Papaloizou 1986; Trilling et al. 1998; 
Trilling, Lunine & Benz 2002). For a disc with surface den- 
sity E(_R, t), the governing equation is, 
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Here, v is the kinematic viscosity which models angular mo- 
mentum transport within the disc gas, and E m is a term 
which allows for mass to be lost from the disc - for exam- 
ple as a consequence of photoevaporation. The second term 
within the brackets describes how the disc responds to the 
planetary torque, A(R, a), where this function is the rate of 
angular momentum transfer per unit mass from the planet 
to the disc. For a planet in a circular orbit at radius a, we 
take, 
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where q = M p /M*, the mass ratio between the planet and 
the star, 



A p = max(H, \R — a\), 



(3) 



and H is the scale height of the disc. Guided by detailed 
protoplanetary disc models (Bell et al. 1997), we adopt H = 
0.057?. 

The transfer of angular momentum leads to orbital mi- 
gration of the planet at a rate, 



(4) 
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if the only torque on the planet comes from the gravitational 
interaction with the disc. 

In the core accretion model for giant planet formation 
(e.g. Pollack et al. 1996), the accretion of the gaseous en- 
velope is predicted to take longer than any other phase of 
the formation process (several Myr in the baseline model 
for Jupiter presented by Pollack et al. 1996). In particular, 
the time scale for accretion is much longer than the time 
scale on which a sufficiently massive planet can establish 
a gap in the protoplanetary disc, since numerical simula- 
tions show that an approximately steady-state gap can be 
set up by a massive planet within ~ 10 2 orbital periods 
(e.g. Lubow, Seibert & Artymowicz 1999). A consequence 
of the inequality of these time scales is that massive planets 
- those of several Jupiter masses - probably accrete most 
of their envelopes subsequent to the development of a gap 
in the protoplanetary disc. Numerical simulations show how 



this accretion may occur. Gas from the outer disc penetrates 
the leaky tidal barrier created by the planet, and flows in- 
ward to form a small circumplanetary disc around the grow- 
ing planet (Lubow, Seibert & Artymowicz 1999; D'Angelo, 
Henning & Kley 2002; Bate et al. 2003). 

The existence of mass flow across gaps onto planets is 
intrinsically a two (or three) dimensional phenomenon (e.g. 
the discussion in Artymowicz & Lubow 1996). The torque 
function (Eq. |5J used in our one dimensional code estab- 
lishes a clean gap for all planet masses above about O.lMj, 
and this gap precludes any mass flow across the gap, or onto 
the planet. To allow for the mass growth of planets, we have 
therefore modified the one dimensional treatment to explic- 
itly include mass flow from the outer disc on to the planet. 
We begin by making an approximate fit to the results of 
two-dimensional numerical simulations (Lubow et al. 1999; 
D'Angelo et al. 2002). We define the efficiency of mass accre- 
tion across the gap via a parameter e, which is the planetary 
accretion rate as a fraction of the disc accretion rate at large 
radii (away from the location of the planet). The results of 
the aforementioned numerical simulations can then be ap- 
proximated by the formula, 
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where Mj is the mass of Jupiter and e ma x is an adjustable 
parameter which can be used to test how the results depend 
upon the overall efficiency of planetary accretion. We use the 
above equation to calculate at each timestep the appropri- 
ate planetary accretion rate. We then remove the required 
amount of mass from the first zone on the outer edge of the 
gap, and add it to the mass of the planet. Note that we as- 
sume that all the mass flow onto the planet originates from 
the outer disc, and do not permit any material to 'bypass' 
the planet and flow directly from the outer disc to the inner 
disc. 

Mass accretion across the gap onto a planet may also be 
expected to lead to accretion of angular momentum, given 
that the specific angular momentum of gas at the outer gap 
edge exceeds that of the planet. This effect - which it is 
easy to show can have a significant influence on the mi- 
gration rate - cannot be straightforwardly measured from 
existing numerical simulations 2 . For this paper, we adopt 
the simplest approach, and assume that the accreted gas 
has the same specific angular momentum as gas at a radius 



Ra 



1.6a. This fixed radius approximates the location of 



the outer edge of the gap throughout most of the calculation. 

Eq. £Q| is solved on a fixed, non-uniform mesh using 
standard explicit numerical methods (e.g. Pringle, Verbunt 
& Wade 1986). The mesh is uniform in a scaled variable 
X oc \fR. Typically, 300 grid points are used, with an inner 
boundary at 0.1 AU and an outer boundary at 200 AU. 
A zero-torque (E = 0) boundary condition is applied at 
Rin- For the protoplanetary disc model adopted, the outer 
boundary is at sufficiently large radius that the choice of 
boundary condition there has no influence on the results. 



2 Bate et al. (2003), for example, explicitly exclude this advcc- 
tive torque from their estimates of the migration rate, due to 
difficulties in measuring the net torque from the small scale cir- 
cumplanetary disc formed in their simulations. 
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2.2 Protoplanetary disc model 

We model the protoplanetary disc as a viscous accretion disc 
(Lynden-Bell & Pringle 1974) suffering mass loss at large ra- 
dius as a consequence of photoevaporation (e.g. Shu, John- 
stone & Hollenbach 1993). Motivation for considering mod- 
els of this form is provided first and foremost by observations 
of photoevaporative flows in Orion (Johnstone, Hollenbach 
& Bally 1998), and is discussed further by Clarke, Gendrin 
& Sotomayor (2001), Matsuyama, Johnstone & Hartmann 
(2003), and Armitage, Clarke & Palla (2003). Briefly, we 
write the viscosity as a fixed (in time) power-law in radius. 
For most of the calculations we have adopted the form, 



v = 3 x 10 1 



(l Au) 



(6) 



This yields a steady-state surface density profile E oc R~ x , 
which is similar to that derived from more detailed proto- 
planetary disc models over the radii of interest (Bell et al. 
1997). To test how sensitive the conclusions are to this as- 
sumption, we have also run one model with a different form 
for the viscosity, 

3/2 



1.3 x 10 1 



(l Au) 



Mass loss from the disc scales with radius as, 



E^i — 
E.,,, OC 



R < Rg 

R > R Q 
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with R g — 5 AU. We express the normalization of the mass 
loss via a parameter M w i nd , which is defined as the total 
mass loss from the disc for a disc with an outer edge at 25 
AU. The instantaneous rate of mass loss will therefore differ 
from this value depending upon the extent of the disc. 

The initial surface density profile for the v oc R disc 
model is, 



E = E 1 




R/R a 



(9) 



while the v oc R?l 2 model is identical except for the re- 
placement of 1/72 by 1/R 3 / 2 . Here, Eo is a constant used 
to define the initial accretion rate, while Ro is a truncation 
radius which sets a smooth exponential cut off to the surface 
density at large radius. For R <C Ro, a disc described by this 
initial condition has a constant accretion rate, so we spec- 
ify the initial surface density of our models via this inner 
accretion rate Minit- 



2.3 Results 

Fig. shows how the evolution of the accretion rate and 
disc mass varies with the strength of photoevaporative mass 
loss. We have computed models with our standard viscos- 
ity (v oc R) that have M wind = 10" 9 M^yr' 1 , M wind = 
2.5 x 10" 9 M©yr _1 , and M wind = 5 x 10" 9 M Q yr _1 . For 
consistency with observational determinations of protoplan- 
etary disc parameters in nearby star-forming regions (e.g. 
Gullbring et al. 1998), we adopt for all of these models an 
initial accretion rate of Minit = 5 x 10 -8 Moyr -1 , and a 
truncation radius of Rq — 10 AU. This yields an initial 
disc mass of 0.066 Mq. An additional model (shown as the 
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Figure 1. Evolution of the accretion rate (upper panel) and mass 
(lower panel) of the protoplanetary disc models used for migra- 
tion calculations. The solid curves show the evolution for models 



10" 



M Q yr- 



with M wind = 10" 9 Mgyr" 1 , M wind = 2.5 X 
and Af wind = 5 X 10 — 9 M0yr _1 (with increasing mass loss rates 
leading to smaller lifetimes). The dashed curve shows a variant 
model with v oc R 3 / 2 and M wind = 5 X 10" 9 Mgyr" 1 . The other 
parameters of the models are as described in the text. 



dashed curve in Fig. was calculated with the v oc 

R 3/2 

viscosity law and a mass loss of M win< j = 5 x 10~ 9 M Q yr _1 . 
With identical choices of Minit and Ro, the initial disc mass 
for this model was 0.067 Mq. 

As expected from previous calculations (Clarke et al. 
2001; Matsuyama et al. 2003), all four models show quali- 
tatively similar evolution. There is an initial phase in which 
the disc mass and accretion rate decline slowly, due primar- 
ily to mass accretion onto the star. Subsequently, the mass 
and accretion rate drop more rapidly as the evolution be- 
comes dominated by mass loss via the wind (Clarke, Gen- 
drin & Sotomayor 2001). Higher rates of mass loss reduce 
the disc lifetime, but all models have observationally accept- 
able lifetimes in the range between 4 Myr and 7 Myr. Most 
importantly for our purposes, in all four models mass loss 
is at least reasonably important (relative to accretion) for 
the overall evolution of the disc. The fraction of the initial 
disc mass that is lost in the wind varies between 40 percent 
and 52 percent for the v oc R models, and is 39 percent for 
the v oc J? 3//2 model. We note that the fraction lost in the 
wind is only a weak function of M w i nd , because the inter- 
val of time over which the wind acts is reduced for higher 
instantaneous mass loss rates. 

To investigate how planets migrate within the evolving 
disc, we run the protoplanetary disc models repeatedly. In 
each run, we add an initially low mass planet (0.1 Mj) into 
the disc at a specified time and radius. We then allow the 
planet to grow and migrate within the evolving disc, and 
record the final planet mass and orbital radius once the disc 
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Figure 2. The final orbital radius and final mass of planets following disc migration, shown as a function of the planets' formation 
epoch. The dashed and solid curves show results for planets formed at initial orbital radii of 5 AU and 10 AU respectively. The left-hand 
panels show the extent of migration in a model disc with M w ; n( j = 10 -9 Mqyt -1 , the centre panels M w ; n( j = 2.5 X 10 — 9 MQyr - x , while 
the right-hand panels depict results from the M w ; nc [ = 5 X 10 — 9 M0yr _1 model. 



has been either accreted or lost in the wind. By varying 
the formation time £f orm , and the formation radius af orm , we 
study how the final outcome depends upon when and where 
in the disc massive planets form. 

Fig. [5] shows the results for the disc with the standard 
viscosity law and varying rates of mass loss. We considered 
planet formation radii of 5 AU and 10 AU, safely outside 
any estimate of the snow line (Sasselov & Lecar 2000), and 
took e max = 1- The final planet mass and orbital radius 
are plotted as a function of the formation time. For the 
low and intermediate rates of mass loss, the sense of orbital 
migration is predominantly inward. Planets formed near the 
end of the disc lifetime end up in orbits close to where they 
formed, accrete relatively little disc gas, and remain as low 
mass objects. Planets formed earlier migrate inwards under 
the action of gravitational torques, and have time to grow 
to several Jupiter masses. These results are consistent with 
previous studies of migration (Trilling et al. 1998; Armitage 
et al. 2002; Trilling, Lunine & Benz 2002). 

For the highest rate of mass loss from the disc, how- 
ever, qualitatively different evolutionary tracks, shown in the 
right-hand panels of Fig. [3] are obtained. The enhanced mass 
loss means that there is a larger range of disc radii across 
which the radial velocity of the gas is outward. This can 
drive significant outward migration. For our choices of pa- 
rameters, we find that inward migration persists for all plan- 
ets formed at 5 AU, while outward migration is the rule for 
planets formed at 10 AU. For formation times between about 



1 Myr and 2 Myr, migration approximately stalls (similar to 
the behaviour reported by Matsuyama, Johnstone & Murray 
2003, though for slightly different reasons), while for earlier 
formation times the gas can drive these outer planets to radii 
of around 20 AU or greater. Significant accretion onto the 
planet occurs throughout this time, so the planets stranded 
at larger radii are all predicted to be massive objects. 

Any attempt to distil the inherently multi-dimensional 
physics of planetary migration into a fast one dimensional 
scheme is bound to be approximate, and there are particu- 
larly obvious uncertainties in our models for planet growth 
and disc evolution. We have already demonstrated, as our 
main result, that for planets forming at radii of around 
10 AU a switch between inward and outward migration oc- 
curs when M w i n d is varied by a factor of a few. The mass loss 
rate via photoevaporation is clearly a vital control parame- 
ter. To check how important some of the other parameters 
are, we have recalculated the migration of planets formed 
at 10 AU in two different models. In one, we altered the 
assumed disc viscosity (to v oc R 3 ^ 2 rather than v oc R), 
with the new viscosity chosen to produce a steady-state sur- 
face density profile E oc R~ 3 ^ 2 . This scaling is one used 
often in studies of the Solar nebula (Weidenschilling 1977). 
A second model was computed with the standard viscosity 
law, but with an accretion efficiency parameter e max = 0.5. 
This change halves the rate of growth of planets via accre- 
tion. Both models used the high rate of mass loss previously 
found to be conducive to outward migration. 
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Figure 3. Sensitivity of the migration results to changes in the 
model parameters. The solid curve shows the extent of migration 
in the standard disc model with v oc R, e max = 1, and Af w ; n( j = 
5 X 10 — 9 MQyr - 1 . The short dashed curve shows a model with 
a different viscosity law (y oc ij 3//2 ), and the long dashed curve a 
model with less efficient accretion on to the planet (e max = 0.5). 



Fig. [3] shows the final radii attained by planets in the 
three models. Outward migration occured in all three mod- 
els, reflecting the primary importance of the assumed mass 
loss rate in determining the fates of the model planets. 
Changing the efficiency of accretion on to the planet made 
negligible difference to the final planetary radii (though it 
reduced the final masses of the model planets by approx- 
imately a factor of two). Substantially greater migration, 
however, was obtained in the calculation with the different 
viscosity law, despite the fact that a smaller fraction of the 
initial disc mass was actually lost to the wind in this case. 
We interpret this as being a side effect of the different surface 
density profile. Outward migration occurs when the torque 
from the inner disc exceeds that from the disc at radii be- 
yond the planet's orbit. The steeper surface density profile 
of the E oc R~ 3 ^ 2 means that there is less mass initially 
exterior to the planet. As this mass is lost in the wind, the 
now unbalanced torque from the inner disc is more effective 
in driving the planet outward. 

2.4 Observational implications 

What are the implications of our gas disc migration calcu- 
lations for the origin of planets at large orbital radii? We 
believe that three general conclusions can be drawn. First, 
mass loss from the outer disc can drive substantial outward 
migration, even when the mass loss is modest enough that 
the disc can survive for several Myr. Higher rates of mass loss 
would lead to more dramatic migration, but the resulting 
short disc lifetimes might preclude planet formation. In the 
central regions of the Orion Nebula, for example, Johnstone 



et al. (1998) infer mass loss rates between 2 x 10 -8 M0yr _1 
and 6 x 10 -7 M0yr _1 , with correspondingly small estimated 
disc lifetimes. Our results suggest that photoevaporation can 
have a major impact on planetary migration even in substan- 
tially more benign environments. Second, outward migration 
driven by the gas disc is favoured in systems where photoe- 
vaporative mass loss is stronger. The predicted rate of mass 
loss due to photoevaporation for a Solar mass star is rather 
small (Shu, Johnstone & Hollenbach 1993), so for relatively 
isolated stars we would expect outward migration only for 
stars significantly more massive than the Sun. An alternative 
possibility is that the mass loss is driven by external irra- 
diation, as in the case of Orion (Johnstone, Hollenbach & 
Bally 1998) . Finally, outward migration via this mechanism 
is a relatively slow process, which occurs on the viscous time 
of the protoplanetary disc. There is ample time for initially 
low mass planets to accrete substantial gaseous envelopes, 
so we would expect planets at large radius to be massive 
objects. 



3 MIGRATION IN MULTIPLE PLANET 
SYSTEMS 

3.1 Introduction 

If the initial outcome of the planet formation process is a 
system of several massive planets, subsequent gravitational 
interactions can lead to many possible outcomes. Planets 
may collide, be ejected from the system, or settle into quasi- 
periodic or periodic orbits. Theoretical investigations of the 
orbital evolution of two-planet systems (Ford, Havlickova, & 
Rasio 2001) and few-planet systems (Chambers, Wetherill, 
& Boss 1996) have been conducted with detailed descrip- 
tions of special cases, such as equal mass planets on coplanar 
orbits. 

Gladman (1993) established that systems with two close 
planets exhibit chaotic but quasi-periodic behavior given the 
appropriate initial conditions. Ford, Havlickova, and Rasio 
(2001) explored the evolution of such systems when the mass 
of each of two planets revolving around a solar-type star in 
nearly circular, coplanar orbits is equal to 10 _3 Mq. We will 
expand on this study by considering different mass ratios for 
the two planets revolving around a central star. In doing so, 
we will show that outward migration is possible, for either 
planet. 



3.2 Methods 

3.2.1 Motivation 

The setup for our set of simulations is motivated primar- 
ily by the gravitational scattering experiments performed 
by Ford, Havlickova, & Rasio (2001). That study consid- 
ered the interaction between two planets of mass 
(on the order of a Jupiter mass) and in an initially close 
configuration. The resulting branching ratios of system out- 
comes was explored. Systems became unstable by ejecting a 
planet or by a collision, one between planets or one between 
a planet and a star. Stable systems remained or settled into 
quasi-periodic orbits over 2 Myr. 

The initial configuration of both planets in Ford, 
Havlickova, & Rasio's (2001) simulations is motivated by 
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Table 1. Classifying the extent of migration or lack thereof for different combinations of mass ratios, with initial parameters < e± < 0.01, 
< e2 < 0.01, 0° < ii < 5°, 0° < «2 < 5°, randomly chosen initial orbital angles, and initial outer semimajor axis that lies within the 
region of quasi-periodic orbits. "Inner" implies that the inner planet exhibited most of the outward migration, "outer" implies that the 
outer planet exhibited most of the outward migration, "both" implies that both planets exhibited outward migration, "none" implies 
neither planet exhibited significant migration, and "*" implies that < 10% of the 300 systems run were stable. 
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Figure 4. Schematic illustration of the initial conditions for scat- 
tering experiments. We randomly populate orbits in the shaded 
region, which lie outside the global chaos boundary but are inside 
the boundary defining guaranteed Hill stability. 



the "Hill Stability Criterion" , an analytic result first applied 
to planetary systems by Gladman (1993). A system is said 
to be Hill stable if the planets cannot approach each other 
closely for all time. We adopt Gladman's (1993) notation: A 
represents the least separation for which both planets, for 
sure, will be Hill stable. Nothing can be said about the Hill 
stability of a system for a separation less than A. A c ^ rep- 
resents the greatest separation at which "global chaos" will 
occur. A system that is "globally chaotic" might produce 
collisions or ejections. For a separation greater than A c ^, 
the planets might exhibit stable quasi-periodic orbits. The 
geometry is illustrated in Fig. [I] Gladman's (1993) analytic 
second order expansion for small values of and fj.2 yields 
an approximate expression for A, 
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Figure 5. Scatter plots of the final eccentricity vs. semimajor axis 
for both planets in the stable systems where m = 5 X 10 — 4 , (j.2 = 
1 X 10" 3 . 



Ford, Havlickova, & Rasio (2001) ran simulations where 
the initial semimajor axis separation of the planets lay be- 
tween A and A c h, so that the planets would neither be in a 
Hill Stable configuration nor become unstable immediately. 
In this work, we sample the entire initial separation range 
spanned by A and A c h in order to best detect planets that 
may migrate outwards and remain on quasi-periodic orbits. 



2 • 3f(/ii +M2) 3 + 2 ■ 33(^1 +/x 2 ) 3 



11/il + 7/J.2 



3~ (pi + ^ 2 ) 3 

Unlike A, an approximation for A c h can only be found 
empirically. By developing overlap resonance criteria for 
two-planet systems, Wisdom (1980) derived the following 
approximation when /ii = /12 = /i: 

A ch ^ 2/J 



3.2.2 Simulation Setup 

We denote the semimajor axis of the initially inner planet 
Oi, the initially outer planet 02, and the respective planet 
mass/central mass ratios as (ii and \i2- Other orbital ele- 
ments will distinguish the planets with a subscript of "1" or 
"2". All runs were performed with ai = 1 AU, so that the 
results may be scaled easily to any ratio of semimajor axes, 
such as for ai = 5 — 10 AU. 

All numerical runs of the three-body system were per- 
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formed with a Bulirsch-Stoer routine from the HNbody inte- 
gration package (Rauch &: Hamilton 2002). The routine ran 
for 2 Myr with an accuracy parameter of 10 -12 , and was re- 
run with uniformly smaller initial timesteps when necessary 
until this accuracy was achieved. For each run, the Bulirsch- 
Stoer routine began with an initial timestep of 0.05 yr, and 
the orbital elements of each body were output every 0.01 
Myr. In most cases, energy and angular momentum errors, 
expressed by (E — Eo)/Eo and (Lz — Lz )/Lz , where Eo 
is the initial total system energy and Lz is the initial to- 
tal angular momentum in the direction perpendicular to the 
orbital plane, did not exceed 10~ 7 . In the most pathologi- 
cal cases, energy and z-angular momentum errors were con- 
served to within 10~ 4 . Angular momenta in the other two 
directions were typically conserved to two order of magni- 
tudes better than the z-angular momentum. 

We are most interested in the evolution of stable sys- 
tems in which both planets remain bound. We define stable 
systems as systems where both planets have semimajor axes 
which never exceed 10 3 AU and eccentricities which never 
exceed 0.99. Further, we deemed a system unstable if at any 
time the HNbody code outputted a negative value for the 
semimajor axis or eccentricity of either planet. 

A preliminary exploration of phase space revealed that 
of the systems that become unstable, most did so within 
0.05 Myr. Thus, each system underwent several checks at 
0.05 Myr; unstable systems were terminated, and stable sys- 
tems were evolved for a total of 2 Myr. In order to reduce 
computer time, we imposed additional conditions on systems 
after 0.05 Myr. Systems not satisfying these conditions were 
terminated. Given that ai = 1 AU for each run, the condi- 
tions are 

1) 1.2 AU < a 2 < 3.0 AU, 

2) ai < 1.05 AU, 

3) If a 2 /ai < 1.5, then e 2 < 0.1, 

4) If a 2 /ai > 1.5, then e 2 < 0.3. 

These conditions were chosen based on a preliminary explo- 
ration of the properties of unstable systems. For example, 
we found that after 0.05 Myr, if the semimajor axis of the 
outer planet is within 50% of oi, then a high (> 0.1) ec- 
centricity of the outer planet implies that the system will 
become unstable. 

By definition, long-term chaotic behavior may differ 
drastically due to an infinitesimal change in initial con- 
ditions. Therefore, because computers use finite-precision 
arithmetic, individual runs are largely irreproducible from 
machine to machine. In this context, one can only make 
statements about the probability of a system behaving in a 
certain manner (Quinlan & Tremaine 1992). Thus, for each 
particular dimension of phase space explored, we performed 
300 runs, each with randomly chosen initial orbital param- 
eters that lay within specific ranges. 

We expanded on the results of Ford, Havlickova, and 
Rasio (2001) by investigating the effects of altering the 
planet /star mass ratio for each planet. For each pair of mass 
ratios, 300 runs were performed, each run with both planets 
having a randomly chosen initial eccentricity between and 
0.01 and initial inclination between 0° and 5°, along with 
an argument of perihelion, longitude of ascending node, and 
mean anomaly each randomly chosen between 0° and 360°. 

We performed 9 sets of 300 runs for each of the following 
four outer planet mass ratios pti: lx 10~ 4 , 5x 10~ 4 , 1 x 10 -3 , 



and 5 x 10 . Each set of runs corresponds to a different 
inner planet mass ratio. We varied the mass of the inner 
planet within one order of magnitude of the mass of the 
outer planet. As mentioned, the initial semimajor axis of 
the inner planet for every run was set at 1.0 AU. The initial 
semimajor axis of the outer planet was chosen in a regime 
that exhibited unpredictable behavior but allowed for quasi- 
periodic orbits. Explanation of this regime follows. 

As mentioned previously, the region bounded by A 
and A c h describe systems that may exhibit Hill instabil- 
ity and may exhibit quasi-periodic motion. When both are 
exhibited, some planets exhibit outward migration, which 
is the topic of this section. So as to include all regions 
from which migrating behavior might arise, we took A c h = 

2 

2 • [min(fii , (12)] 7 ■ Because ai = 1 AU for each run, values 
of a 2 were randomly chosen between 1 + A ch and 1 + A for 
the simulations. 

Our initial conditions for the scattering experiments are 
similar to those adopted by most previous authors (e.g. Ford, 
Havlickova & Rasio 2001; Marzari & Weidenschilling 2002; 
Chambers, Wetherill & Boss 1996), in that we start with 
fully-formed planets in relatively close proximity to each 
other, and consider only their mutual gravitational inter- 
actions. Such initial conditions are chosen primarily for sim- 
plicity and compatibility with previous work, and can only 
partially be justified on physical grounds. In particular, since 
many of the trial systems prove to be unstable over time 
scales comparable to the lifetime of protoplanetary discs 
(Haisch, Lada & Lada 2001), there is likely to be an ear- 
lier epoch during which both planet-planet gravitational in- 
teractions and planet-disc interactions are important. This 
earlier stage of evolution could lead to an unstable multi- 
ple planet system in at least two ways. First, two planets 
might form in well-separated orbits, but then subsequently 
migrate into an unstable configuration as a consequence of 
disc-driven migration (Goldreich & Tremaine 1980; Lin & 
Papaloizou 1986; Lin, Bodenheimer & Richardson 1996). For 
this to happen, the planets would have to avoid becoming 
trapped into resonance during the migration process (Snell- 
grove, Papaloizou & Nelson 2001; Lee & Peale 2002; Mur- 
ray, Paskowitz & Holman 2002). Alternatively, the planets 
might form close together, and be stabilized against imme- 
diate violent instability by the presence of a surrounding 
gas disc (Lin & Ida 1997; Nagasawa, Lin & Ida 2003) 3 . Fur- 
ther, unstable planetary systems may form after the disc has 
dissipated. Planet-planet interactions alone in crowded sys- 
tems of Jovian-mass planets most often leaves only two sur- 
vivors in close, quasiperiodic but ultimately unstable orbits 
(Adams & Laughlin 2003) . Although any of these pathways 
could plausibly lead to initial conditions that are similar to 
those which we (and other authors) have assumed, a full 
treatment will obviously need to model the more complex 
interactions that are possible during the phase when the gas 
disc is being dispersed. 



3 In related work, Agnor & Ward (2002) considered the damping 
of terrestrial planet eccentricity by a remnant gas disc. 
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Figure 6. Cumulative probability distributions for stable systems with (12 = 1 x 10 -4 . For clarity, values of fi\ are abbreviated such 
that 1 X 10~ 5 = 1(— 5). The upper panel represents orbital parameter distributions of the inner planet, and the lower panel the outer 
planet. 



3.3 Results 

Table 1 and Figs. 10191 displays the data we took. Table 1 
summarizes the data in terms of mass ratio and extent of 
migration, while the figures provide a more detailed look at 
migration behavior. Fig. provides a representative sample 
of the final semimajor axes and eccentricities of both planets 
when significant outward migration occurs, and provides a 
perspective on how data is plotted in Figs. 10191 For exam- 
ple, the 43 stable systems displayed in Fig.|^|are plotted as 



the solid lines in Fig. |S] Fig. |5] illustrates that planets that 
migrate both inward and outward tend to increase their ec- 
centricity. Although the extent of migration in both direc- 
tions is similar, those planets that migrate inward tend to 
increase their eccentricities at a higher rate. 

Figs. l619l each represent a set of six cumulative probabil- 
ity distributions that display the fraction of stable systems 
vs. semimajor axis ratios, eccentricities, and inclinations of 
the outer and inner planet after exactly 2 Myr. The upper 
panel in each figure represents the initially outer planet, and 
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Figure 7. Cumulative probability distributions for stable systems with (12 = 5 X 10 -4 . For clarity, values of fi\ are abbreviated such 
that 1 X 10~ 4 = 1(— 4). The upper panel represents orbital parameter distributions of the inner planet, and the lower panel the outer 
planet. 



the lower panel represent the initially inner planet. Each fig- 
ure keeps the mass of the outer planet fixed, but varies the 
mass of the inner planet. Although we obtained nine curves 
for each graph on each figure, we show only the curves that 
exhibit significant radially outward migration, plus some 
that do not, and we do not show any curves where less than 
10% of the 300 runs were stable. 

Detecting migration for the outer planets is difficult 
because their semimajor axes were randomly chosen for 
each run. To aid in determining the net radial movement 



of these planets, we define a.\ = a± (final)/ai (initial) and 
02 = a-i (final) /a? (initial). Each panel in Figs. lol§l contain 
distributions of semimajor axis ratios, rather than absolute 
semimajor axes, because those are scalable in this study. 

Table 1 and Figs. IbliJl illustrate that outward migration 
of both the inner and outer planets occurs more frequently 
with smaller values of fi2- Fig. |S| displays extensive migra- 
tion. In the lowest mass case, almost 60% of the stable inner 
planets migrate at outward to at least 150% of their original 
semimajor axis. As ji\ is increased, migration occurs less fre- 
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Figure 8. Cumulative probability distributions for stable systems with (12 = 1 x 10 -3 . For clarity, values of fi\ are abbreviated such 
that 5 X 10~ 4 = 5(— 4). The upper panel represents orbital parameter distributions of the inner planet, and the lower panel the outer 
planet. 



quently; for fii = 2.5 x 10~ 4 no inner planets migrated out 
to 150% of their original semimajor axis. The correlation of 
mass ratio to final state behavior is less apparent but still 
present in the final eccentricity and inclination curves. In 
both graphs, the least-mass case prompts the highest final 
eccentricity and inclination values for the orbit of the inner 
planet. 

The final state of the outer planet, represented by the 
lower panel of Fig. |SJ reflects the exchange of angular mo- 
mentum in this system. One sees that when an inner planet 



migrates outward, it leaves the initially outward planet back 
so that Qf2 < 1. The one instance where the outer planet ex- 
periences significant migration is the highest mass case - the 
same case where cvi < 1.5. The eccentricity and inclination 
distributions for the outer planet in this massive case dif- 
fer drastically from the behavior seen for other inner planet 
masses. 

In Fig. |7| /12 = 5 x 10 -4 and the result is significantly 
less migration than seen in Fig. © Similarly to Fig. © the 
highest mass case exhibits the greatest migration of the 
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Figure 9. Cumulative probability distributions for stable systems with (12 = 5 X 10 -3 . For clarity, values of fi\ are abbreviated such 
that 7.5 X 10 — 3 = 7.5(— 3). The upper panel represents orbital parameter distributions of the inner planet, and the lower panel the outer 
planet. 



outer planet. In this case, ~20% of the outer planets and 
~15% of the inner planets migrate beyond 150% of their 
original semimajor axis. Halving the inner mass value re- 
duces the probability for ai < 1.5 or Q2 < 1.5 to ~5%. In 
contrast to Fig. the least mass case fails to show signifi- 
cant outward migration of the inner planet, despite the inner 
planet exhibiting the highest final inclinations. 

In Fig. |S] where the outer mass is a full order of magni- 
tude more massive than in Fig.|S| we see even less migration. 
Inner bodies for the smallest mass case exhibit high eccen- 



tricities and inclinations, but only ~5% satisfy ct\ > 1.5. 
The case where fii — /xa = 1 X 10 -3 is the same case studied 
by Ford, Havlickova, & Rasio (2001). Their Figs. 11-12 can 
be compared to the semimajor axis ratio distributions, ec- 
centricity distributions, and inclination distributions in Fig. 

m 

No migration of the inner planet occurs for any case in 
Fig. |5] However, the outer planet shows significant migra- 
tion. Further, as shown by the sequence of four curves in 
the lower panel, the higher the mass of the outer planet, 
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Figure 10. Cumulative probability distributions for a system with negligible migration (/ti = (12 = 5 X 10 — 3 , upper panel) and for a 
system with significant migration (fix = 1 X 10~ 5 , (12 = 1 x 10 — 4 , lower panel) for different tolerances, which are given with abbreviated 
scientific notation. 900 runs were performed for the upper panel system, and 300 were performed for the lower panel system. Similarly 
to Figs. ISI91 only the stable systems are shown. 



the greater the extent of the migration. For the higher 
mass cases, about 25% of the initally outer planets satisfy 
Q2 > 2.5. Remarkably, these planets all retain final eccen- 
tricities < 0.4, in stark contrast to the large final eccentric- 
ities of the migrating planets of Figs. 10181 

Because of the chaotic nature of the three body prob- 
lem, one encounters difficulty when predicting the appropri- 
ate timescale over which to integrate. We adopted a simi- 
lar timescale to the one used by Ford, Havlickova, & Rasio 
(2001), but recognize that systems that appear to be stable 
at 2 Myr might become unstable at some future time. We 
extended the running time of one system to sample the con- 
sequences. After 2 Myr, 55 out of the 300 systems for the 
case fjbi = 5 X 10 -3 , [12 = 1 X 10~ 3 remain stable. By run- 
ning this system for 10 Myr, we found only 3 out of the 55 
systems became unstable, and did so quickly. Further, these 
three systems became unstable before 3 Myr. 

For any set of 300 runs of a chaotic system, the slightest 
change in any input parameter of the Bulirsch-Stoer algo- 
rithm might drastically alter the results of any individual 
runs, but keep the same global behavior. The extent of the 
invariance of this global behavior is a function of the code 
used and the number of systems sampled. To explore this 



measure of invariance for the runs presented here, we reran 
two systems with four different initial tolerances (5 x 10 -11 , 
1 x 10~ n , 5 x 10~ 12 , 1 x 10~ 12 ). One of these systems 
(fii — (12 — 5 x 10 -3 ) exhibited negligible migration, and 
the other (fix = 1 X 10~ 5 , /12 = 1 x 10~ 4 ) exhibited sig- 
nificant migration. In order to achieve the largest feasible 
sample size for this error analysis, we tripled the number 
of runs performed to 900 for the system exhibiting negli- 
gible migration. As Fig. 1101 illustrates, runs with different 
tolerances are practically indistinguishable for the case of 
no migration, but vary up to 20% for the case of significant 
migration. 



3.4 Summary 

The gravitational interaction of a pair of unequal mass plan- 
ets that lie in close initial configurations allow either planet 
to migrate outward to at least twice it's initial semima- 
jor axis. Planets that migrate outward lie in quasi-stable, 
slightly inclined eccentric orbits with an eccentricity that 
spans the entire permissible range and an inclination up to 
about 10°. Although either planet in any given system may 
drift outward, the smaller the mass of the planets, the higher 
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the tendency for the initially inner planet to migrate out- 
ward. Further, less massive giant planets such as Uranus or 
Neptune are more likely to migrate outward than Jupiter- 
mass planets. This result, derived from our gravitational 
scattering simulations alone, is consistent with Thommes, 
Duncan & Levison's (1999, 2002) conclusion that Uranus 
and Neptune's current location is a result of outward migra- 
tion amongst Jupiter and Saturn. 



4 CONCLUSIONS 

As a first step in explaining the presence of planets at large 
orbital radii, we have shown that outward migration of pro- 
toplanets is possible both by planet-disc interactions and by 
planet-planet gravitational scattering without the presence 
of a disc. Strong mass loss in discs coupled with planetary 
cores that are formed at about ~ 10 AU allow planets to mi- 
grate outward in discs to radii that are as much as a factor 
of several in excess of their initial semimajor axes. Planets 
that migrate in such a manner are likely to be massive. We 
predict that gas-driven outward migration should be most 
likely to occur around more massive stars, whose strong UV 
flux can drive a powerful photoevaporative outflow. Sub- 
sequently, when in the appropriate chaotic regime, planets 
within a multiple planet system may migrate outward due to 
gravitational scattering alone. Planets that migrate in this 
manner may be massive or not, however low mass objects 
tend to exhibit the most extensive outward migration. Or- 
bital migration is typically accompanied by an increase in 
eccentricity that spans the allowable range for elliptic orbits. 
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